2020 excess mortality & voting patterns in CH
Modelling redistributed community deaths
Data
Spatial
Downscaled data
Prepared in 08.Rmd.
exp_deaths_2020_year = read_rds("results/exp_deaths_2020_year.Rds") %>%
select(-cant_exp_deaths, -cant_observed, -p) %>%
mutate(id_kt = as.integer(as.factor(canton))) Removing communities without voting data.
exp_deaths_2020_year = exp_deaths_2020_year %>%
filter(!is.na(vote_yes_nov_cat)) %>%
filter(!is.na(vote_yes_jun_cat)) Removing <40 age group.
Create special variables for inla.
exp_deaths_2020_year = exp_deaths_2020_year %>%
mutate(id_space=as.numeric(as.factor(GMDNR)),
id_space2=id_space,
sex_fem=ifelse(sex=="Female",1,0),
age_num=as.numeric(as.factor(age_group)),
age_60s=ifelse(age_group=="60-69",1,0),
age_70s=ifelse(age_group=="70-79",1,0),
age_80s=ifelse(age_group=="80+",1,0),
type_urban=ifelse(r_urban1=="Urban",1,0),
type_rural=ifelse(r_urban1=="Rural",1,0),
density_high=ifelse(r_urban2=="Dense",1,0),
density_low=ifelse(r_urban2=="Low",1,0),
sep5=ifelse(median_ssep3_q=="5th - highest",1,0),
sep4=ifelse(median_ssep3_q=="4th",1,0),
sep3=ifelse(median_ssep3_q=="3rd quintile",1,0),
sep2=ifelse(median_ssep3_q=="2nd",1,0),
sep1=ifelse(median_ssep3_q=="1st - lowest",1,0),
lang_fr=ifelse(r_lang=="French",1,0),
lang_it=ifelse(r_lang=="Italian",1,0),
vote_nov_q5=ifelse(vote_yes_nov_cat==5,1,0),
vote_nov_q4=ifelse(vote_yes_nov_cat==4,1,0),
vote_nov_q3=ifelse(vote_yes_nov_cat==3,1,0),
vote_nov_q2=ifelse(vote_yes_nov_cat==2,1,0),
vote_nov_q1=ifelse(vote_yes_nov_cat==1,1,0),
vote_jun_q5=ifelse(vote_yes_jun_cat==5,1,0),
vote_jun_q4=ifelse(vote_yes_jun_cat==4,1,0),
vote_jun_q3=ifelse(vote_yes_jun_cat==3,1,0),
vote_jun_q2=ifelse(vote_yes_jun_cat==2,1,0),
vote_jun_q1=ifelse(vote_yes_jun_cat==1,1,0),
E=log(ifelse(munici_exp_deaths==0,1e-4,munici_exp_deaths)))Modelling observed and expected
Step 1: model development
To facilitate model development we only use the median excess by municipality, age group and sex in 2020.
data1 = exp_deaths_2020_year %>%
group_by(canton, GMDNR, GMDNAME, age_group, id_space, sex, munici_observed, munici_pop,
density_high, density_low, across(starts_with("sep")), border, lang_fr, lang_it,
across(starts_with("vote")),type_urban,type_rural) %>%
summarise(munici_exp_deaths=median(munici_exp_deaths),
munici_excess=median(munici_excess)) %>%
mutate(E=log(ifelse(munici_exp_deaths==0,1e-4,munici_exp_deaths)))
rm(exp_deaths_2020_year)
gc() used (Mb) gc trigger (Mb) max used (Mb)
Ncells 4342919 232.0 7852384 419.4 7852384 419.4
Vcells 7603824 58.1 1036210404 7905.7 1292298099 9859.5INLA setup
Model 1.0: no covariates
We use a model structure similar to Poisson regression, where \(O_{t,i,j,k}\), the number of observed deaths during week \(t\) in municipality \(i\), age group \(j\) and sex group \(k\), depends on the number of expected deaths \(E_{t,i,j,k}\) based on historical data and a linear predictor \(\lambda\).
\[ O_i \sim \text{Poisson}(E_i \times \exp(\lambda)) \\ \] At start, the linear predictor \(\lambda\) only includes one intercept parameter \(\alpha\), so that the estimate of \(\exp(\alpha)\) can be interpreted as an average relative excess mortality for 2020. By adding covariates to \(\lambda\), we aim to disentangle the various factors that are associated with excess mortality at the local level.
We implement this model in R-INLA, a Bayesian inference package that is especially adapted to spatial data. This is achieved in practice by including \(\log (E_{i,j,k})\) as an offset (although an alternative formulation based on the E argument exists). During model development, we compare different model versions based on the WAIC (lower values imply a better fit).
model1.0 = INLA::inla(munici_observed ~ 1 + offset(E),
data = data1,
family = "Poisson",
control.compute = list(config = TRUE, waic = TRUE),
quantiles = c(0.025, 0.5, 0.975),
num.threads = threads,
safe = TRUE)
summary(model1.0)
Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts,
", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
verbose, ", " lincomb = lincomb, selection = selection, control.compute
= control.compute, ", " control.predictor = control.predictor,
control.family = control.family, ", " control.inla = control.inla,
control.fixed = control.fixed, ", " control.mode = control.mode,
control.expert = control.expert, ", " control.hazard = control.hazard,
control.lincomb = control.lincomb, ", " control.update =
control.update, control.lp.scale = control.lp.scale, ", "
control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
num.threads, ", " blas.num.threads = blas.num.threads, keep = keep,
working.directory = working.directory, ", " silent = silent, inla.mode
= inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
.parent.frame)")
Time used:
Pre = 0.589, Running = 0.754, Post = 0.11, Total = 1.45
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 0.325 0.004 0.317 0.325 0.332 0.325 0
Watanabe-Akaike information criterion (WAIC) ...: 49278.96
Effective number of parameters .................: 3.53
Marginal log-Likelihood: -24648.51
is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)') mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 1.383534 1.003654 1.373679 1.383534 1.39346 1.383534 1[1] 1.383562As a sanity check, we find a relative excess mortality of 38% for 2020, that is coherent with a simple calculation (74,776 observed / 54,046 expected = 1.38). Remember that we excluded the age group 0-40, which explains why this is higher than numbers reported for Switzerland, generally around 10% for 2020. We can also look at the model fit and at the residuals. Obviously the model fit is not good here, as this basic model assumes a unique relative excess mortality for all areas, sexes and age groups.
Model 1.1: age and sex
We hypothesize that excess mortality affected different age and sex groups differently. We thus add the age group, the sex and the interaction of the two as covariates.
model1.1 = INLA::inla(munici_observed ~ - 1 + offset(E) +
sex:age_group,
data = data1,
family = "Poisson",
control.compute = list(config = TRUE, waic = TRUE),
quantiles = c(0.025, 0.5, 0.975),
num.threads = threads,
safe = TRUE)
summary(model1.1)
Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts,
", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
verbose, ", " lincomb = lincomb, selection = selection, control.compute
= control.compute, ", " control.predictor = control.predictor,
control.family = control.family, ", " control.inla = control.inla,
control.fixed = control.fixed, ", " control.mode = control.mode,
control.expert = control.expert, ", " control.hazard = control.hazard,
control.lincomb = control.lincomb, ", " control.update =
control.update, control.lp.scale = control.lp.scale, ", "
control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
num.threads, ", " blas.num.threads = blas.num.threads, keep = keep,
working.directory = working.directory, ", " silent = silent, inla.mode
= inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
.parent.frame)")
Time used:
Pre = 0.654, Running = 0.996, Post = 0.112, Total = 1.76
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
sexFemale:age_group40-59 0.008 0.024 -0.040 0.008 0.055 0.008 0
sexMale:age_group40-59 0.354 0.018 0.318 0.354 0.390 0.354 0
sexFemale:age_group60-69 0.035 0.020 -0.003 0.035 0.074 0.035 0
sexMale:age_group60-69 0.417 0.015 0.388 0.417 0.447 0.417 0
sexFemale:age_group70-79 0.235 0.013 0.210 0.235 0.260 0.235 0
sexMale:age_group70-79 0.473 0.011 0.453 0.473 0.494 0.473 0
sexFemale:age_group80+ 0.493 0.006 0.481 0.493 0.504 0.493 0
sexMale:age_group80+ 0.149 0.007 0.136 0.149 0.163 0.149 0
Watanabe-Akaike information criterion (WAIC) ...: 47190.67
Effective number of parameters .................: 4.17
Marginal log-Likelihood: -23651.68
is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)') mean sd 0.025quant 0.5quant 0.975quant
sexFemale:age_group40-59 1.007864 1.024454 0.9612512 1.007864 1.056738
sexMale:age_group40-59 1.424950 1.018531 1.3745805 1.424950 1.477166
sexFemale:age_group60-69 1.036069 1.019763 0.9970809 1.036069 1.076581
sexMale:age_group60-69 1.517741 1.015056 1.4739324 1.517741 1.562852
sexFemale:age_group70-79 1.264505 1.012778 1.2334258 1.264505 1.296368
sexMale:age_group70-79 1.605326 1.010613 1.5724496 1.605326 1.638889
sexFemale:age_group80+ 1.636631 1.006044 1.6174162 1.636631 1.656075
sexMale:age_group80+ 1.161138 1.007045 1.1452712 1.161138 1.177224
mode kld
sexFemale:age_group40-59 1.007864 1
sexMale:age_group40-59 1.424950 1
sexFemale:age_group60-69 1.036069 1
sexMale:age_group60-69 1.517741 1
sexFemale:age_group70-79 1.264505 1
sexMale:age_group70-79 1.605326 1
sexFemale:age_group80+ 1.636631 1
sexMale:age_group80+ 1.161138 1[1] -2088.29As expected, the relative excess mortality varies a lot across age and sex groups. It’s very small in females aged 40-59 and 60-69 (in fact the data is compatible with no excess in both cases). It increases in females aged 70-79, and even more so aged 80+. It’s comparatively higher in males below 80, but somewhat surprisingly lower in males in age group 80+. We observe an improvement of the model fit, not easy to spot on the plot because of the large number of points, but made clear by the large decrease in WAIC.
Model 1.2: spatial variability
We now account for spatial variability, first in a simple way using an i.i.d. random effect, so that all municipalities can vary independently from each other around a global average.
model1.2 = INLA::inla(munici_observed ~ - 1 + offset(E) +
sex:age_group +
f(id_space, model = "iid"),
data = data1,
family = "Poisson",
control.compute = list(config = TRUE, waic = TRUE),
quantiles = c(0.025, 0.5, 0.975),
num.threads = threads,
safe = TRUE)
summary(model1.2)
Call:
c("inla.core(formula = formula, family = family, contrasts = contrasts,
", " data = data, quantiles = quantiles, E = E, offset = offset, ", "
scale = scale, weights = weights, Ntrials = Ntrials, strata = strata,
", " lp.scale = lp.scale, link.covariates = link.covariates, verbose =
verbose, ", " lincomb = lincomb, selection = selection, control.compute
= control.compute, ", " control.predictor = control.predictor,
control.family = control.family, ", " control.inla = control.inla,
control.fixed = control.fixed, ", " control.mode = control.mode,
control.expert = control.expert, ", " control.hazard = control.hazard,
control.lincomb = control.lincomb, ", " control.update =
control.update, control.lp.scale = control.lp.scale, ", "
control.pardiso = control.pardiso, only.hyperparam = only.hyperparam,
", " inla.call = inla.call, inla.arg = inla.arg, num.threads =
num.threads, ", " blas.num.threads = blas.num.threads, keep = keep,
working.directory = working.directory, ", " silent = silent, inla.mode
= inla.mode, safe = FALSE, debug = debug, ", " .parent.frame =
.parent.frame)")
Time used:
Pre = 0.61, Running = 2.1, Post = 0.238, Total = 2.95
Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
sexFemale:age_group40-59 0.011 0.024 -0.036 0.011 0.058 0.011 0
sexMale:age_group40-59 0.358 0.018 0.322 0.358 0.394 0.358 0
sexFemale:age_group60-69 0.039 0.020 0.001 0.039 0.078 0.039 0
sexMale:age_group60-69 0.420 0.015 0.391 0.420 0.450 0.420 0
sexFemale:age_group70-79 0.238 0.013 0.213 0.238 0.263 0.238 0
sexMale:age_group70-79 0.477 0.011 0.456 0.477 0.498 0.477 0
sexFemale:age_group80+ 0.497 0.006 0.485 0.497 0.509 0.497 0
sexMale:age_group80+ 0.153 0.007 0.139 0.153 0.167 0.153 0
Random effects:
Name Model
id_space IID model
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for id_space 3453.45 2510.49 1339.89 2741.78 9707.17 2116.51
Watanabe-Akaike information criterion (WAIC) ...: 47159.19
Effective number of parameters .................: 14.49
Marginal log-Likelihood: -23646.50
is computed
Posterior summaries for the linear predictor and the fitted values are computed
(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)') mean sd 0.025quant 0.5quant 0.975quant
sexFemale:age_group40-59 1.011058 1.024506 0.9641761 1.011058 1.060222
sexMale:age_group40-59 1.430421 1.018607 1.3796297 1.430420 1.483087
sexFemale:age_group60-69 1.039892 1.019829 1.0006107 1.039891 1.080718
sexMale:age_group60-69 1.522570 1.015131 1.4783840 1.522569 1.568084
sexFemale:age_group70-79 1.268978 1.012874 1.2375430 1.268977 1.301220
sexMale:age_group70-79 1.610709 1.010717 1.5773914 1.610707 1.644743
sexFemale:age_group80+ 1.643565 1.006297 1.6234816 1.643556 1.663948
sexMale:age_group80+ 1.165188 1.007210 1.1488960 1.165184 1.181732
mode kld
sexFemale:age_group40-59 1.011057 1
sexMale:age_group40-59 1.430419 1
sexFemale:age_group60-69 1.039890 1
sexMale:age_group60-69 1.522567 1
sexFemale:age_group70-79 1.268974 1
sexMale:age_group70-79 1.610703 1
sexFemale:age_group80+ 1.643538 1
sexMale:age_group80+ 1.165177 1[1] -31.47579The age and sex effect remains similar, but the model fit as measured by the WAIC is improved now that we account for local differences. We can observe this municipality effect, that applies in all age and sex groups of the municipality in exactly the same way.